#ifndef __METHODSNUM_H
#define __METHODSNUM_H 1

//#include <feel/options.hpp>
//#include <feel/feelcore/parameter.hpp>
//#include <feel/feelcore/feel.hpp>
//#include <feel/feelalg/backend.hpp>

//#include <functionSup.cpp>

//#include <research/life-workspace/Fluid-Structure/cldata.hpp>
#include <research/life-workspace/Fluid-Structure/cltype.hpp>

#include <research/life-workspace/Fluid-Structure/functionSup.cpp>



namespace Feel {


    template </*typename cl_type,*/typename Appli>
    class MethodsNum
    {
    public :

        typedef Appli appli_ptrtype;
        typedef typename appli_ptrtype::element_type appli_type;

        typedef MethodsNum<appli_ptrtype> self_type;

        typedef typename appli_type::backend_ptrtype backend_ptrtype;

        typedef typename appli_type::vector_ptrtype vector_ptrtype;
        typedef typename appli_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;


        typedef typename appli_type::element_dirichlet_type element_dirichlet_type;

        appli_ptrtype M_appli;

        backend_ptrtype M_backend;

        MethodsNum(Appli __app)
            :
            M_appli(__app),
            M_backend(appli_type::backend_type::build( M_appli->application()->vm()))
        {}


        template <typename cl_type>
        void
        linearSolver(cl_type & __cl, vector_ptrtype U)
        {
            auto condlim = __cl;
            auto Xh = M_appli->functionSpace();

            vector_ptrtype F( M_backend->newVector( Xh ) );
            sparse_matrix_ptrtype A( M_backend->newMatrix(Xh,Xh) );

            M_appli->updateLinearPDE(U,A,F);

            std::cout << "[MethodsNum] : linearSolver solve start\n";
            M_backend->solve( _matrix=A,
                              _solution=U,
                              _rhs=F,
                              _maxit=1000
                              //_atolerance=1e-50,
                              //_dtolerance=1e-15
                              );
            std::cout << "[MethodsNum] : linearSolver solve finish\n";

        }

        template <typename cl_type>
        void
        AlgoNewton(cl_type & __cl,vector_ptrtype U)
        {

            std::cout << "[MethodsNum] : Newton Algo start\n";

            auto condlim = __cl;
            auto Xh = M_appli->functionSpace();
            auto dirElem = M_appli->getDirichletElement();

            ////////////////////////////////////////////////////////////////////////

            if (true)//!M_appli->weakCL())
                {
#if 0
                    ForEachCL( condlim,cl::dirichlet_vec,
                               modifVec222(*dirElem,U, Expression , PhysicalName ) );
                    ForEachCL( condlim,cl::paroi_mobile,
                               modifVec222(*dirElem,U, vf::idv(M_appli->meshVelocity()) , PhysicalName ) );
#else
                    M_appli->updateCLDirichlet(U);
#endif
                }
            ////////////////////////////////////////////////////////////////////////
            vector_ptrtype Uold( M_backend->newVector( Xh ) );
            vector_ptrtype R( M_backend->newVector( Xh ) );
            sparse_matrix_ptrtype J( M_backend->newMatrix(Xh,Xh) );
            vector_ptrtype DeltaU( M_backend->newVector( Xh ) );

            double nUold;
            double Err=1;
            double ErrTol=1e-8;//8;//1e-5;
            uint cpt=0;

            M_appli->updateResidual( U, R );
            while (Err>ErrTol && cpt<1000)
                {
                    Uold->zero();
                    Uold->add(1.,*U);
                    nUold=Uold->l2Norm();

                    R->scale(-1.0);

                    M_appli->updateJacobian( U, J, R );


                    std::cout << "[MethodsNum] : Newton solve start\n";
                    M_backend->solve( _matrix=J,
                                      _solution=DeltaU,
                                      _rhs=R,
                                      _maxit=10000,
                                      _atolerance=1e-50,
                                      _dtolerance=1e7
                                      );
                    std::cout << "[MethodsNum] : Newton solve finish\n";
                    *U += *DeltaU;
                    //M_appli->updateCLDirichlet(U);
                    M_appli->updateResidual( U, R );

#if 1
                    Err=R->l2Norm();
#else
                    Uold->add(-1.0,*U);
                    if (nUold!=0) Uold->scale(1./nUold);
                    if (cpt>0) Err=Uold->l2Norm();
#endif

                    //std::cout<<"[MethodsNum] : Newton Residu L2 : " << Err << "\n";
                    std::cout<<"[MethodsNum] : Newton  : iter " << cpt << " Residu L2 : " << Err << "\n";

                    /*
                     *(M_appli->getSolution()) = *U;
                       M_appli->exportResults();
                    */

                    ++cpt;

                }
            //M_appli->updateCLDirichlet(U);
            std::cout << "[MethodsNum] : Newton Algo finish\n";
        }


        template <typename cl_type>
        void
        AlgoNewton2(cl_type & __cl,vector_ptrtype U)
        {

#if 1
            auto condlim = __cl;
            auto Xh = M_appli->functionSpace();
            auto dirElem = M_appli->getDirichletElement();

            ////////////////////////////////////////////////////////////////////////

            ForEachCL( condlim,cl::dirichlet_vec,
                       modifVec222(*dirElem,U, Expression , PhysicalName ) );
            ForEachCL( condlim,cl::paroi_mobile,
                       modifVec222(*dirElem,U, vf::idv(M_appli->meshVelocity()) , PhysicalName ) );

            ////////////////////////////////////////////////////////////////////////

            vector_ptrtype R( M_backend->newVector( Xh ) );
            sparse_matrix_ptrtype J( M_backend->newMatrix(Xh,Xh) );
            vector_ptrtype DeltaU( M_backend->newVector( Xh ) );

            //M_backend->nlSolve( J, Uvec, R, 1e-10, 10 );
            //boost::shared_ptr<SolverNonLinear<double> >
            //    nlsolver( SolverNonLinear<double>::build( SOLVERS_PETSC ) );

            //M_backend->nlSolver()->setType(TRUST_REGION);
            // set up the non linear solver
            /*
            M_backend->nlSolver()->residual = boost::bind( &NavierStokes<Dim, Order, Entity>::updateResidual,
                                                           boost::ref( *this ), _1, _2 );
            M_backend->nlSolver()->jacobian = boost::bind( &NavierStokes<Dim, Order, Entity>::updateJacobian,
                                                           boost::ref( *this ), _1, _2, _R );
            */
            M_backend->nlSolver()->residual = boost::bind( &appli_type::updateResidual,
                                                           boost::ref( *M_appli ), _1, _2 );
            M_backend->nlSolver()->jacobian = boost::bind( &appli_type::updateJacobian,
                                                           boost::ref( *M_appli ), _1, _2, R );

            M_backend->nlSolver()->setRelativeResidualTol(1e-16);
            M_backend->nlSolver()->setAbsoluteResidualTol(1e-16);
            M_backend->nlSolver()->setNbItMax(100);
            M_backend->nlSolver()->solve( J, U, R, 1e-10, 100 );
#endif

        }


    };


} // end namespace Feel


#endif //__METHODSNUM_H
